DETERMINATION OF ^O NT ACT PARAMETERS 
USED IN DEM MODELS FOR REALISTIC 
SIMULATION OF BALL MILLS 


A Thesis Submitted 
in Partial Fulfilment of the Requirements 
for the Degree of 
MASTER OF TECHNOLOGY 


by 

ASUTOSH GOPINATH 



DEPARTMENT OF MATERIALS & METALLURGICAL ENGINEERING 
INDIAN INSTITUTE OF TECHNOLOGY KANPUR 
APRIL, 1999 



1 8 MAY 1999 

CENTRAL ubRart 

1, >. T., KANPUR 


I lllllll III! mil Hill lll'l 11,11 llli 1,11 


CERTIFICATE 



It is certified that the work contained in this thesis titled 11 Determination of 
Contact Parameters used in DEM Models for Realistic Simulation of Ball Mills", 
by Asutosh Gopinath, has been carried out under my supervision and that this 
work has not been submitted elsewhere for a degree. 


Date: Dr - 

(Associate Professor), 

Materials and Metallurgical Engineering, 
Indian Institute of Technology, Kanpur, 
April, 1999. 



B. K. 


Dedicated to.. 


my parents 



Acknowledgements 


As my academic career at I. I. T. Kanpur is coming to an end with the submis- 
sion of this thesis, I take the opportunity to express my sincere gratitude to my 
respected guide Dr. B. K. Mishra for the enlightenment he has bestowed on me 
about the interesting field of computer simulation. The successful completion of 
this work has been possible only due to his excellent guidance throughout this 
work. 

I would like to thank all my friends for providing me a cordial and friendly 
environment and thus making my stay at I.I.T. Kanpur a memorable one. 



Abstract 


The discrete element method (DEM) has come a long way in simulating multi- 
ple interacting deformable bodies undergoing large absolute or relative motions. 
These separate bodies can also undergo progressive fracturing and result in the 
generation of smaller distinct bodies. Therefore, in order to correctly predict the 
motion and fracture behavior of these distinct elements it is essential to accu- 
rately determine the forces developed at the contact points. In the DEM the 
contact forces are determined by describing the contact behavior through a pair of 
spring and dashpot. The responses of the spring and the dashpot to any applied 
external force are the subject of the present investigation so as to make them 
realistic. 

DEM has been used sucessfully to simulate the motion of balls in tumbling 
mills. While the overall motion of the balls appear to be realistic, a careful look 
at the corresponding quantitative features such as power draw and impact energy 
distribution reveal that there is a scope to improve upon the contact model. What 
is needed is a unified framework within which the contact parameters would be 
specified such that for any given system these parameters remain insensitive to 
any change within the system. 

In the present work, at first, a bouncing ball problem is considered in detail. 
The motion of the ball is described by a non-linear differential equation that 
takes into account a typical spring-dashpot type contact behavior. This equation 
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is solved numerically to completely describe the motion of the ball before and 
after the impact. The corresponding experimental work reported in the literature 
on a very sensitive ultra-fast load cell provided a basis for correlation. Thus, it 
has become possible to directly extract the contact parameters that are to be 
embedded in the analytical model. These parameters include the contact stiffness 
reflecting material and geometric properties of the two contacting surfaces, and 
the damping coefficient which is a measure of the coefficient of restitution. 

In the DEM simulation, identifying the critical time-step is very crucial, as 
a larger time-step can lead to unrealistic or erroneous results, and a small time- 
step may give more accurate results but requires large computer memory and 
a long simulation time. It was found that the contact parameters obtained by 
non-linear model necessitated a small critical time-step. One way of overcoming 
this large computational effort was examined using the "Equivalent Linearisation" 
technique. By using this method an equivalent linear model was obtained from 
the non-linear model and equivalent parameters were calculated. Interestingly, the 
equivalent contact parameters so obtained required a larger time-step than in the 
non-linear model, which significantly reduced the overall computational effort. 

A comparative investigation of the linear, the non-linear and the equivalent 
linearized models was done for the bouncing ball problem. It was found that 
linearized model was more accurate than the linear model and required lesser 
computational effort than the non-linear model. 

In the concluding part of the present work, DEM. simulation of a laboratory 
size tumbling mill is performed to compare the power draw obtained for the three 
models. 
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Chapter 1 


Introduction and Literature Review 


Discrete element method (DEM) is a computational technique to simulate the 
behavior of any system comprising discrete interacting particles. In this method, 
the contact forces acting between particles are represented by a pair of spring and 
dashpot. For any two or more interacting particles within the system, one pair 
of spring and dashpot is activated for each contact when the distance between 
the centre of particles become less than sum of their radii. The contact forces 
and other body forces are summed to determine the net force acting on a given 
particle. Then the net acceleration of the particle is determined by using Newton’s 
second law of motion. The positions of the particles are then determined using an 
explicit numerical time integration recipe such as 


&t+At = it + x t At 

(1.1) 

Xt+At = x t + (x t + x t+ A t) 0.5 At 

(1.2) 


Thus, this numerical scheme allows simulation of discretely interacting bodies in 
an explicit manner that has a wide range of engineering applications. 

The DEM that describes the motion of assemblies of particles was pioneered by 
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Cundall and Strack [1], Cundall and Strack proposed a model consisting of spring, 
dashpot and slider components at contact points of adjacent particles. This model 
is commonly used for calculating the particle impact force which is then used for 
updating the particle position. Several assumptions about the spring and damping 
forces used in impact calculations have been made by different researchers. Until 
few years back linear spring and dashpot model was used to represent a contact, 
but lately it has been found that a nonlinear spring and dashpot model gives 
more accurate results [2]. The accuracy of simulation also depends on the accurate 
determination of contact parameters which facilitate a better estimate of the inter- 
elemental forces. 

The representation of contact as a pair of spring and dashpot has enabled DEM 
to explore those areas in which interactions between discrete particles of various 
size and shapes is taking place. Over the past two decades extensive research on 
DEM has been carried out and it has been successfully applied to simulate discrete 
particle interactions in many areas. Table 1.1 [3] outlines the areas of application 
of discrete element method. 

In many DEM applications in general, and in ball mill in particular, it has been 
found that the material parameters selected to simulate spring and dashpot play 
the most critical role. Even a slight variation in the value of material parameters 
is likely to give a large variation in the simulation results. Therefore, when DEM 
is used to predict the power draft of a ball mill which is considered to be a good 
performance indicator, it is imperative that the parameters are determined care- 
fully. These parameters essentially embody the properties of the system: stiffness, 
damping, and friction. The contact is typically represented by a system of spring 
and dashpot as shown in Figure 1.1. The forces in the spring and the dashpot 
as a result of the contact is determined by the associated stiffness and damping 
quantities. In the following few paragraphs a brief description on work done by 



Table 1.1: Areas of application of the discrete element method. 

Application Areas 

•Mineral processing 

- Simulation of ball-charge motion in tumbling mills 

- Simulation of bed stratification in jigs 

• Behavior of soils, rocks and granular materials 

Rock pit slope stability behavior 

Underground structure stability in jointed rocks 

Study of earthquake mechanism and plate tectonics 

Modeling of geological structure evolution and formation of reservoir traps 
Liquification under dynamic loading 
Blast survivability studies 
Hazardous waste disposal studies 

• Impact and explosive dynamics 

Automobile crash simulation 

Projectile penetration and explosive containment 

Kinetic energy weapon effects 

Vulnerability and particle tracking studies during explosive fragmentation 

• Mechanical behavior 

Metal forming 

Interaction of machinery components 
Analysis of linkages and chains 
Vibration control and feedback studies 
Fracture mechanics 
Study of joints and limbs 
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Spring 

Slider 


Dashpot 


Figure 1.1: Particle contact model (disk-disk contact). 


various researchers in modeling the contact and. the various ways of determining 
the contact parameters is presented. 

There are basically four types of force-displacements relationships that are 
used to describe the contact behavior [2]. These relationships and the associated 
methods for determining the contact parameters are described below: 


• Linear contact law is the simplest one in which the contact parameters are 
independent of the history of contact forces and are typically represented as 


spring 


= kx 


(1.3) 


Fdashpot — Q% (f-4) 

where, k is the stiffness parameter and q is the viscous damping coefficient. 
The spring and dashpot are simply placed between the contacting particles to 
simulate the load-displacement responses. Since the contact parameters are 
assumed to be constant, this law fails to simulate the pressure-dependency 
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of parameters of the medium and thus is not realistic. 

• Non-linear contact law is the second method for determining contact param- 
eters. In this method, a non-linear representation of contact force is done. 
The contact parameters are represented as 

E spring — (F^) 


F iashpot — F 


( 1 . 6 ) 


where r is the non-linearity in spring and s is the non-linearity in dashpot. 
The contact parameters are determined by using Mindlin solution for con- 
tact. For the case of identical elastic rough spheres of radius r s and elastic 
constants E s and v s , the normal stiffness k can be written as 


_ ' 2 E S A 

(1 - *.) 


( 1 . 7 ) 


where A is the radius of contact area. The tangential stiffness, on the other 
hand, depends upon the complete history of normal and shear forces. This 
law gives a more realistic simulation but requires large computational time 
and memory. 

• Modified linear contact law was proposed to save computer time and mem- 
ory, and to avoid the complication of Mindlin solution. In this method, the 
normal stiffness was taken same as in Mindlin solution, but the shear con- 
tact stiffness was made a linear function of current normal force acting at the 
contact. This aided in doing away with the cumbersome task of determine 
the complete history of normal and shear forces which might vary from sys- 




Section 1.1 Contact Force Analysis for Spring Dashpot 
Models 


6 


tern to system. The main limitation of the law is for cyclic load simulation, 
where it does not predict the realistic behavior of contact. 

•- Modified non-linear contact law was proposed to reduce the computational 
time and memory requirement. In this method, the normal stiffness is the 
Mindlin stiffness give by Eq. 1.7, and the shear stiffness is some non-linear 
function of the current value of normal and tangential contact forces. 

The four contact laws outlined above serve as a basis for modeling contacts in 
various discrete particle systems. Quite a few models have been developed either 
directly using the above laws or by slightly modifying them. The following sections 
deals with work done in modeling system comprising discrete particles using the 
above mentioned contact laws. 


1.1 Contact Force Analysis for Spring Dashpot 
Models 

The linear contact model was first used by Cundall and Strack [1]. This is the 
most commonly used model that assumes a linear spring and a linear damper. 
The force is applied when the particles overlap, that is when the distance between 
their centres is less than the sum of the particle radii. For a linear spring and 
dashpot, the equation of motion of an impacting particle can be expressed as: 

ma + qct + ka = 0 (1.8) 

where m is the mass of ball, a is the particle overlap, q is the damping coefficient, 
and k is the contact stiffness. 

In order to determine contact stiffness, Zhang et al [4] limited the maximum 
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anticipated inter-particle penetration, <y max , to be a small fraction of the particle 
diameter d. This constraint was achieved by defining: 

k = f 2 mv 2 0 /d 2 (1.9) 


where v Q is the anticipated maximum speed of any particle in the system, and / 
is the penetration factor defined by putting a max equal to dj f In this study the 
value of / was chosen to be equal to 10.0. The damping coefficient q was defined 
by: 


q = 2 In (1/e) 


kmiirij/ ( mi + mj ) 
7r 2 + (In (1/e))' 


1 !/2 


( 1 . 10 ) 


where k is the contact stiffness, and rrii and rrij are the masses of the disks i and 
j respectively. 

For the discrete particle simulation of two-dimensional fluidized bed, Tsuji et al 
[5] used the linear model by assuming a constant value for stiffness. Qualitatively, 
the results were satisfactory in many respects, such as particle circulating motion 
and mixing, but the model could not predict the correct force versus time trend 
(Figure 1.2). 

Experiments carried out by Zhang and Whiten [6] have shown that the expres- 
sion for damping is incorrect, and added that the damping coefficient q should 
not be constant, but is expected to be a function of displacement. They showed 
(Figurel.3) that when the damping is non-zero, the initial force is large, which is 
contrary to experimental measurements that show the force rising smoothly from 
zero to maximum value and then returning smoothly to zero. 

Over the last few years, quite a few papers have been published by different 
researchers reporting modeling of contact force by non-linear spring and dashpot. 



Normalised force m Pressure (Pa) 
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Time (sec) 


: Pressure fluctuation: (a) experiment; (b) calculation [5]. 


Normalised time 


0 0.5 1 1.5 2 2.5 3 3.5 4.0 



Figure 1.3: Normalized force verses normalized time using linear model, a is nor- 
malized damping [6]. 
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The non-linear impact phenomenon was first and extensively studied by Hertz. For 
linearly elastic and homogeneous, with perfectly smooth contact surfaces, Hertz 
[7] calculated that the force acting at the contact area (known as Hertz spring 
force) using: 


F SV ring = ka 3/2 ( 1 . 11 ) 

where k is the material stiffness and a is the elastic approach of the colliding 
bodies. Hence, the equation of motion of head-on colliding bodies can be written 
as: 


ma + ka 3 ^ 2 = 0 


( 1 . 12 ) 


where, m is the mass of body. 

It was experimentally proved by Velusami [7] that a contact damping force 
also existed at the impact area. According to him, the damping force is given 
by Fdamping — qoi 3 l 2 a so, the equation of motion of the impacting body can be 
written as: 


ma -1- qa 3 ^ 2 a + ka 3 ^ 2 = 0 


(1.13) 


where, q is the damping constant. 

Hertz theory is limited to metal-metal impact only. Realizing the implication 
of Hertz’s analysis, researchers in the field of DEM applied non-linear force theory 
to obtain a more realistic model for contact force. Yet another non-linear model 
came up, but this time by Tsuji et al [8]. But the model used to represent damping 
differed greatly with that of Velusami. Here the the damping coefficient was found 
heuristically and was taken as ( = a(mk) 1/2 a l/4 . The equation of motion was 
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written as 


ma + (a + ka = 0 


(1.14) 


or, 


ma + a{mk) l ^ 2 a l ^a + ka 3 ^ 2 = 0 


(1.15) 


where, a is a constant which depends on coefficient of restitution e. 

The contact stiffness k was determined by contact theory. In the case of two 
spheres of same size (radius equal to r s ), k is expressed by: 


_ \/2r s E s 
~3 JT^uf) 


(1.16) 


where, E s is the Young’s modulus, and v s is the Poisson ratio of the particles. In 
case of contact between a sphere and a wall, k is expressed by: 


k = 


Wf7 

3 

kzd -l lzi& 
E s ' E w 


(1.17) 


where, E w is the Young’s modulus, and u w is the Poisson ratio of the wall. Tsuji 
et al used the above non-linear model to simulate the plug flow of cohesionless 
particles in a horizontal pipe [8]. This model however showed marked quantitative 
disagreement with experiment, but was satisfactory in the sense that the method 
needed fewer empirical factors than the existing analyses of plug flow. 

Goldsmith [9] came up with another non-linear model using the theory of 
compressive strain wave propagation in solid, that described the motion of a high- 
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speed impact of an elastic sphere on the end of a uniform bar as: 

ma + ^ ^- a l/2 a + ka 3 ^ 2 — mg = 0 (1-18) 

2 pAc 0 

where, g is the acceleration due to gravity, p is the density of bar, A is the cross- 
section area of bar, c 0 — is the velocity of propagation of strain wave down 
the bar, and Y is the elastic modulus of bar. 

King and Bourgeois [10] used this model to simulate the energy required for 
particle breakage. The diameter of rod was 3/4 inch, the diameter and mass of 
ball used were 1 inch and 0.067 kg, respectively. The expression for k used here 
was same as that used by Tsuji et al [8] for ball-wall impact. Again, in this case, 
quantitative disagreement regarding variation of contact force versus time was 
observed ( Figure 1.4). 



Figure 1.4: Comparison of predicted and measured responses from 5 cm drop 
height [10]. 
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1.2 Limitations of Earlier Works 

For simulation of a ball mill, a detailed study of the parameters involved in the 
process of grinding is necessary. The most important parameters from the point 
of view of power requirement during grinding are stiffness and damping. All the 
models outlined in this chapter are purely theoretical and are heavily based on 
certain assumptions that might drastically affect their performance .when they are 
directly applied to simulate the behavior of ball-charge motion in industrial mills. 

The assumptions which these models make are given below: 

• In all the above models, the damping term is assumed to be a constant in- 
dependent of the amount of deformation undergone by the particles. This 
results in erroneous calculation of contact forces [6]. Initially, just prior to 
contact the damping force should be zero and then it should change grad- 
ually. Tsuji’s model, on the other hand, is able to do this qualitatively but 
quantitatively it requires some more adjustment. Also, Tsuji’s and Gold- 
smith’s model is inefficient in predicting force-displacement history which is 
very important because it gives an idea about energy loss of the system. 

• Since the DEM deals with individual particles with normal and shear contact 
springs and dashpots, it is not possible to specify overall aggregate properties 
as the Young’s modulus, coefficient of friction, coefficient of restitution, and 
Poisson ratio. Instead, it is necessary to choose approximate values of con- 
tact parameters [12]. The quantitative discrepancies in the results obtained 
by various researchers can be attributed to this. 
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1.3 Scope of Present Work 

In the present work, instead of using the contact parameters that were obtained 
theoretically by researchers, an attempt has been made, for the first time, to ex- 
tract the relevant parameters by correlating the model developed with the exper- 
imentally obtained data. This makes the present work more realistic in the sense 
that a model is developed by analysis of experimental data obtained from literature 
[13, 14]. By doing this certain assumptions pertaining to material characteristics, 
e.g. smooth surface and no prior deformation are done away with. Secondly, the 
damping term is modeled in such a way that the unrealistically high initial force 
is suppressed and quantitative agreement with experiment is achieved. Finally, 
average values of contact parameters which are obtained by using the model de- 
veloped in present work, for different mill conditions, are selected for simulation of 
power draw of a laboratory size ball mill. This ensures incorporation of individual 
contact characteristics to gross ball-charge behavior. 



Chapter 2 


Determination of Contact Stiffness 
and Damping 


In this chapter, a vast amount of experimental data on single ball impact are 
analyzed. The force- displacement as well as the force-time relationship during a 
collision is considered. These experimental data are fitted to a non-linear spring- 
dashpot model of the following form 

m& + qa s a + ka r = 0 (2.1) 

where k is the stiffness parameter which determines the maximum force attained 
during impact, r is the nonlinearity in the spring, q is the damping parameter 
which determines the amount of energy expended during the impact, and s is a 
measure of the non-linearity in the dashpot. Different values for the parameters 
were assumed and the best values were selected that matched the experimental 
data. The data were obtained by conducting drop-ball tests [13] on ultra fast load 
cell device. 
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2.1 Drop-Ball Test 

In order to determine the most realistic values of the contact parameters, an elab- 
orate analysis of numerous experimental data gathered by means of the ultrafast 
load cell equipment [13] ( UFLC ) was carried out. The description of UFLC is 
given in the following section. 


2.2 Description of Ultrafast Load Cell Device 

Ultrafast load cell (UFLC) device is a device to measure impact responses of 
materials under dynamic loading conditions. The loading and fracture of particles 
in tumbling mills happens within a very short period of time. With UFLC, the 
grinding action can be studied in detail by measuring the impact of a ball dropping 
into a bed of mineral particles. The UFLC has a resolution of 10 -5 second. 

In this apparatus, a drop-weight impacts from a certain height on a particulate 
agglomerate situated on top of a steel rod. The sphere release mechanism can be 
adjusted vertically to vary the drop height of the impacting sphere. The UFLC 
utilizes the propagation of elastic waves. The impact of drop-weight on particles 
results in a pressure on top of the rod which causes a compressive strain wave 
to propagate through the rod. This strain wave is detected by solid state strain 
gauges which are glued to the rod. The signal from the strain gauges is amplified 
and stored in a digital storage oscilloscope which is interfaced to a computer. From 
the recorded wave, the force acting on the particle is computed. A schematic of 
this equipment is drawn in Figure 2.1. 
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Figure 2.1: Schematic outline of the ultra-fast load cell 
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2.3 Experimental Data on UFLC 

A series of drop-ball tests were carried out by Hofier [13] and Mishra [14] at the 
Comminution Center of the University of Utah. Two types of data were gathered 
from drop-ball tests: first, impacts in which there was no particle layer laid on the 
rod, i.e., the ball was directly dropped on the rod, and second, impacts in which 
there were one or more layers of particles laid on the rod. One type of data is useful 
because it represents the metallic collisions of ball-ball and ball-wall type in a ball 
mill, and the other type of data is useful for the determination of parameters to 
simulate material breakage in ball mill. For metal-metal impact with no particle 
layers in between, three balls of mass 0.096, 0.252 and 0.647 kg were dropped 
from a height of 0.3 m. The force-time and force-deformation history on the top 
surface of ultra fast load cell rod was detected by the solid state strain gauges 
and recorded in the oscilloscope. Whereas, the other type of data was obtained 
by conducting experiments for different drop-heights, number of particle layers, 
masses of the impacting ball and types of material placed between the impacting 
bodies. In this case the maximum force developed at the contact and energy loss 
due to impact were recorded. The experimental data is presented in Appendix A 
in Table A.l through Table A. 4. 

In this chapter, a methodology is developed to extract the contact parameters 
from the data obtained by means of drop ball tests. 

2.4 Determination of Contact Parameters for Metal- 
Metal Impact 

The contact parameters were determined by solving the model equation (Eqn. 
2.1) in two steps. First, the range of nonlinear parameters r and s were fixed. The 
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nonlinearity in the spring r was kept in the range 1 — 2. Similarly, the nonlinearity 
in dashpot s was kept in the range 0 — 1. Second, approximate values of contact 
parameters k and q were taken as initial guess values. Then the model equation 
was solved using Fourth Order Runge-Kutta method. The force-displacement and 
force-time relationship fitted with the corresponding experimental results. In case 
of any mismatch, the values of all the model parameters were modified and the 
procedure was repeated all over again. The results showing force-time and force- 
deformation history both for the model and experiment in case of metal-metal 
impacts are shown in Figure 2.2 through Figure 2.4. 

It was found that the predicted result matched the experimental data in all 
three cases for r equal to 1.6 and s equal to 0.8. However, the numerical values 
of k and q varied in each case. Table 2.1 summerises the experimental conditions 
and the corresponding numerical results. The values of stiffness k and damping 
q obtained numerically were used to calculate the energy loss due to metal-metal 
impact. The energy loss was obtained by calculating the area enclosed by force- 
deformation loop. The energy of impact obtained by using k and q that are 
extracted from the non-linear model were in good agreement with the measured 
energy for different ball masses. Thus, for the metal-metal impact the model took 
the following form: 


ma + qa°' 8 a -I- fra 1 ' 6 = 0 (2.2) 

This analysis however, points to the fact that a linear contact model is not 
appropriate inasmuch as it fails to predict the force-time Figure 1.3 history. How- 
ever, with proper tuning of the parameters, the impact energy can be matched for 
a given collision. Although non-linear models can predict the force-time history 
during a collision as seen here, the associated parameters turn out to be numer- 
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Figure 2.2: Comparison between experiment [14] and model for ball-flat impact 
for ball of mass 0.096 kg. 
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Figure 2.3:. Comparison between experiment [14] and model for ball-flat impact 
for ball of mass 0.647 kg. 

Table 2.1: k and q values and comparisons between model and experimental [14] 
energy loss for metal-metal impact. 


Experiment No. 

1 

2 

3 

Ball mass (kg) 

0.096 

0.647 

0.252 

Non-linear stiffness k 

4.8el0 

5.2el0 

5.0el0 

Non-linear damping q 




Energy loss model (J) 


0.985 

SEfll 

Energy loss Experimental (J) 

0.20 

1.08 

0.47 

Absolute error (%) 

6.5 

8.8 

3.4 
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Time(sec) 



Figure 2.4: Comparison between experiment [14] and model for ball-flat impact 
for ball of mass 0.252 kg. 
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ically large. For example, as seen from Table 2.1, the non-linear stiffness values 
are on an average of the order of 10 10 . This turns out to be the most significant 
drawback of employing a non-linear contact model in the discrete element method. 
This aspect is discussed in section 2.5 below. 

2.5 Constraint of Critical Time-step 

Consider a body of mass m suspended from a rigid support by a spring of stiffness 
k. The mass is made to oscillate simple harmonically. The time period At (or 
critical time-step) of its oscillation can be shown to be equal to 

At = 2 y/mjk, (2.3) 

In order to accurately estimate the motion of particles using DEM, where each 
contact is represented by a pair of oscillating spring and dashpot, the time-step 
chosen is 1/10 smaller than the critical time-step. This is done to accurately 
capture the motion of the impacting balls. Hence, in problems involving simulation 
by discrete element method, the critical time-step is given by 


At = 0.2 \/rnJk (2.4) 

where m is the mass of the smallest ball. As it is clear that for very large value 
of k obtained by non-linear model, the critical time-step becomes very small. 
Therefore, such a large value of k would require calculations of order 10 10 for each 
ball in 1 second of real time. Therefore, it would require huge memory space and 
very long simulation time. 

The computational constraint of large memory space and time requirements 
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relieved if the value of k is of a lower order. It has been observed that for the 
linear model, the computational time requirement is quite less, but at the same 
time it does not correctly predict force-time and force-deformation histories. This 
discrepancy of the linear model propagates into the DEM algorithm resulting in 
erroneous calculation of energy loss. Hence, in the present work, a linear model 
was obtained from the non-linear one in such a way that the error in calculation 
of energy loss was minimized. In other words, it can be said that the linear model 
obtained was equivalent to the non-linear model in terms of energy loss prediction. 
The technique which is employed to get the equivalent values of contact parameters 
from the non-linear model by converting it into equivalent linearized form is known 
as M Equivalent Linearisation” and is described in the following section. 

2.6 Equivalent Linearisation 

Consider the non-linear equation of the form: 

x + ef(x,x) = F(t,T) (2.5) 

where e is the perturbation parameter and T is the time period. This equation 
can be approximated by: 


X + PeqX + K eq X = F(t, T) 
and the error E in the approximation is: 


( 2 . 6 ) 


E(x, x) = [ef(x, x) - K eq X - /3 eq x 


(2.7) 
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The average of the square of the error over one cycle is: 


(- E 2 (X,X)) T = — J [ef(x,i) ~ K eq X - PeqXfdt 


Minimizing this quantity with respect to the coefficients K eq and P eq , 


° {E% = ^J\ Edt = o 


{E\ = ~ j xEdt = 0 


( 2 . 10 ) 


Therefore, 


rji 

J x[ef(x, x) - K eq X - Peg x]dt = 0 (2-11) 

Jo 


[ exf(x,x)dt — [ K eq x 2 dt- f p eq xxdt = 0 (2.12) 

Jo Jo Jo 


Since, / 0 T xxdt = 0 for a periodic solution, 


f^eq r T 


/q T exf(x,x)dt 

L T x2dt 


(2.13) 


Similarly, 


Peg = 


Io T x 2dt 


(2.14) 


Having minimized the error, it is then neglected. Hence, the governing non-linear 
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equation reduces to an equivalent linear equation: 

x + p eq x + K eq x = F(t, T ) (2.15) 

Note that both fi eq and n eq depend on the solution x(t). Also, the signs of { E 2 ), 

( E 2 ) and ( E 2 ) indicate that by setting the quantities ( E 2 ) and 

(E 2 ) zero, E 2 is minimized. The equivalent linearisation method is most 
suited for systems with small non-linearities. However, it is applied for ball- 
ball and ball-wall contact also to explore its applicability, even though the non- 
linearities in these cases are much more severe. 

Thus, the non-linear model 

ma + qa 08 a + ka 1-6 = 0 (2.16) 

on equivalent linearisation is transformed into 


ma + q eq a + k eq a = 0 


(2.17) 


where, q eq is the linearized viscous damping coefficient, and k eg is the linearized 
stiffness term. Hence, 


and, 


k eq — 


/p T oif(oi,ai)dt 

So a 2 dt 


Qeq — 


/ 0 T af(a,a)dt 


fo “ 2 * 


(2.18) 


(2.19) 


Thus, for metal-metal impact, the different models investigated in this study are 
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presented below 


Linear model 

: mdi + qu n a + k lin Oi = 0 

(2.20) 

Non — linear model 

: ma + qa 0 ' 8 a + to 1 ' 6 = 0 

(2.21) 

Equivalent — linearised model 

: ma + q eq a + k eq a = 0 

(2.22) 


Table 2.2 shows the equivalent linearized values of contact parameters k eq and 
q eq found by the method outlined above and also the predicted energy loss for 
the case of metal-metal impact. The table also shows the stiffness and damping 
parameters taken for linear model. In case of linear model the value of stiffness 
kun was taken as initial slope of experimental force-displacement curves in Figure 
2.2 through Figure 2.4. The damping constant in this case was determined by [16] 
(for derivation see Appendix B) 

(2.23) 

V In e + 7r 2 

where, m is the ball mass and e is the coefficient of restitution. The value of 
coefficient of restitution was taken as 0.8. 

As seen from Table 2.2, the equivalent linearized parameters give more accurate 
prediction than linear parameters. It is also observed that the equivalent value of 
stiffness entails calculation of the order of 10 s which is significantly lower than the 
order 10 10 of stiffness obtained by non-linear model (see Table 2.1). Hence, it can 
be concluded that equivalent linearized model is a good compromise between linear 
and non-linear models as far as accuracy and computational time requirement is 
concerned. 

Furthermore, in tumbling mills grinding particles, metal-metal impacts are 
less likely to occur. More than 90 % of impacts occur in which particle layers are 
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Table 2.2: k eq and q eq values and comparisons between linear model, non-linear 
model, equivalent linearized model and experimental energy losses for metal-metal 
impact. 


Experiment No. 

1 

2 

3 

Ball mass (Kg) 

0.096 

0.647 

0.252 

Linear stiffness ku n 

1.4e8 

1.26e8 

1.35e8 

Linear damping qn n 

519.5 

1279.4 

826.5 

Equivalent stiffness k eq 

1.17e8 

1.72e8 

1.13e8 

Equivalent damping q eq 

984.74 

1912.86 

1361.34 

Energy loss linear model (J) 

0.181 

0.948 1 

0.446 

Energy loss non-linear model (J) 

0.187 

0.985 

0.454 

Energy loss linearized model (J) 

0.185 

0.967 

0.451 

Energy loss Experimental (J) 

0.20 

1.08 

0.47 

Absolute error % linear model 

9.5 

12.2 

5.1 

Absolute error % non-linear model 

6.5 

8.8 

3.4 

Absolute error % linearized model 

7.5 

10.5 

4.04 


trapped between the balls. Therefore, the contact parameters obtained in this 
chapter cannot simulate the actual mill condition because they represent metal- 
metal impact only. Hence, more analysis was done to obtain the parameters 
in order to simulate the actual mill conditions. The model adopted to analyze 
the force-time and force-displacement behavior during metal-metal impact was 
assumed to hold good where the layers of particles are caught in between the 
impacting bodies. This assumption was necessary because the experimental force- 
time and force-displacement plots were not available for impacts in which particle 
layers present between the impacting bodies. This has direct application in ball 
mill simulation using DEM. Ball mills are used for grinding particles. The tum- 
bling action of the balls impart the required force on the particles. DEM models 
are used to simulate not only the tumbling actions of the entire ball mass but 
also the resulting distribution of forces and power draw. It is essential that the 
parameters in the DEM model be chosen so as to minimize the computational 
time and accurately predict the power draw. 
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In order to meet the above objectives it is proposed to represent the impact 
where layers of particles are caught between the colliding bodies. In other words, 
it is assumed that a hypothetical layer of particles surrounds the impacting balls 
as shown in the following schematic (Figure 2.5). This idea is explored to obtain 
the contact parameters as discussed in Chapter 3. 



Particle layers caught between 
impacting balls 



Equivalent representation 

1 

Non-linear 

model 

1 

Extract k and q 


Figure 2.5: A schematic showing hypothetical layer of particles surrounding the 
impacting balls. 




Chapter 3 

Determination of Contact 
Parameters for Mill Conditions 


Determination of contact parameters under real milling conditions is conceived in 
a different manner. The earlier approach as described in Chapter 2 is only valid 
under hypothetical situation where a ball mill is operated without any particles in 
it. In real mills, particles are ground by the tumbling action of the ball. Therefore, 
the contact parameters are not exactly the one that would correspond to a metal- 
metal collision. In this chapter, pertinent experimental data from the literature is 
collected. These are drop ball test data with layers of particles present between 
the colliding masses. As before, a non-linear model is fitted to the experimental 
data and the associated parameters of the model are extracted. 

3.1 Analysis of Experimental Data 

The contact parameters obtained for metal-metal impacts cannot be used for tum- 
bling mills because more than 90 % impacts capture layers of particles. Therefore, 
in order to obtain the contact parameters it is essential to collect data under con- 
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ditions that would apply to a ball mill. Holler [13] did drop ball tests in the UFLC 
using layers of particles. However, his data (Table A. 2 through Table A.4) basi- 
cally show variations in maximum impact force developed and energy required in 
breaking particles caught between the balls as a function of drop height, number 
of particle layers, and material type. 

3.1.1 Effect of Drop Height on Energy Loss and 
Impact Force 

It is observed that both the impact force and energy loss increase with drop height. 
The Figure 3.1 and Figure 3.2 depict this effect for particle 5 and 10 layers of 
particles of quartz. It is also seen that the number of layer has negligible effect 
on energy loss but it has noticeable effect on maximum force developed on the 
particles. 



Drop height (m) 

Figure 3.1: Effect of drop height on energy loss due to impact of ball of mass 
0.647 kg on layers of quartz [13]. 
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Figure 3.2: Effect of drop height on impact force due to impact of ball of mass 
0.647 kg on layers of quartz [13]. 

3.1.2 Effect of Material Type on Energy Loss and Impact 
Force 

It is observed that energy loss increases from quartz to limestone and further to 
magnetite. The maximum impact force is also found to increase in the same order. 
Both types of behavior can be attributed to the micro-structural features such as 
presence of pores, cracks and other flaws in materials. Materials having lower 
number of flaw per unit volume have higher strength and consequently require 
higher force to deform and more energy to fracture. Similarly, materials having 
greater flaw density are brittle and require lower force to deform and lesser energy 
to fracture (Figure 3.3 and Figure 3.4). 
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Material type 


Figure 3.3: Effect of material type on energy loss for ball of mass 0.252 kg dropped 
from a height of 0.5 m on 2 layers of particles [13]. 
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Figure 3.4: Effect of material type on maximum impact force for ball mass 0.252 kg 
dropped from a height of 0.5 m on 2 layers of particles [13]. 
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3.2 Determination of Contact Parameters for 
Impacts with Particle Layers in Between 

Various conditions that have an effect on contact parameters are [15]: ball mass, 
impact velocity, material to be ground, and number of particle layers. The data 
gathered from literature depicting the affect of the aforesaid parameters are in 
the present study used as a basis for contact parameters determination. The non- 
linear model obtained for the case of metal-metal impact was solved numerically 
to simulate particle breakage as well. Here, the same procedure as in the case of 
metal-metal impact, was followed to solve the non-linear equation. Initial guess 
values of k and q were taken and the model equation (Eqn 2.1) was solved using 
Fourth Order Runge-Kutta method. The values of the parameters were chosen 
in such a way that the maximum force and energy loss obtained matched with 
the experiments. In case of any mismatch, the values of k and q were fine-tuned 
and the procedure was repeated. The numerical solution provided the contact 
parameters for the non-linear model for various mill conditions mentioned above. 
The model was then linearized to get the equivalent linear contact parameters 
also. 

In the following sections, the values of contact parameters obtained are re- 
ported and how the above mentioned conditions effect the contact parameters is 
discussed. 
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3.3 Effect of Impact Velocity or Drop Height of 
Impacting Ball 


It was found that the value of stiffness increased with drop height. The damping 
term also showed the same trend (see Table 3.1). 


Table 3.1: Effect of drop height on contact parameters for ball of mass 0.647 kg 
and 5 layers of particles. 


Material type 

Drop-height (m) 

0.25 

0.50 

1.0 

Limestone 

k 

0.8e8 

10e8 

15e8 

q 

0.5e6 

1.8e6 

1.78e6 

Ic 

™eq 

2.2e6 

2.1e7 

3.0e7 

Qeq 

2740 

7088 

12165 

Quartz 

k 

2.1e8 

14.1e8 

32.8e8 

q 

0.33e6 

1.23e6 

1.30e6 

keq 

4.7e6 

3.5e7 

6.4e7 

Qeq 

1380 

5411 

6258 


It has been shown by Venkatraman and Narayanan [11] that the impact force 
increases with impact velocity raised to the power 6/5 for ball-ball impact. There- 
fore, the impact force is expected to increase with drop height. Since, stiffness 
bears a correspondence with impact force, it also increases with increase in drop 
height. 

It was observed that collision energy increased with velocity of impact. Gold- 
smith [9j showed that coefficient of restitution decreases with increase in impact 
velocity (see Appendix B) and consequently with increase in drop height. Damp- 
ing coefficient being measure of energy loss during collision, increases with decrease 
in coefficient of restitution and therefore with increase in drop height. Moreover, 
damping coefficient is also a function of stiffness and increases with increase in 
stiffness (see Figure 3.5). 




Equivalent damping Equivalent stiffness 
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Figure 3.5: Effect of drop height on stiffness and damping values of materials for 
ball of mass 0.647 kg and 5 layers of particles. 
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3.4 Effect of Number of Particle Layers 


It was found that the number of particle layers decreased the contact stiffness and 
increased the damping coefficient (see Table 3.2). 

Table 3.2: Effect of layers of particles on contact parameters for ball of mass 
0.647 Kg dropped from the height of 0.25 m. 


Material type 


Layers 

2 

5 

10 

Limestone 

k 

1.6e8 

1.5e8 

1.4e8 

9 

1.2e6 

1.23e6 

1.78e6 

^ eq 

3.9e6 

3.6e6 

3.2e6 

i 

Qeq 

1949 

2740 

3310 

Quartz 

k 

2.38e8 

2.1e8 

2.0e8 

9 

0.5e6 

0.52e6 

0.6e6 

keq 

4.5e6 

4.2e6 

3.5e6 

Qeq 

1000 

1380 

1500 


This is because particles in additional layers increase the cushioning effect on 
the contact points and reduce the amount of maximum force developed and hence 
reduce the stiffness. 

Also the energy absorbing capacity of layers increase due to the presence of 
more particles between the balls, thus increasing damping (see Figure 3.6). 


3.5 Effect of Ball Mass 

The ball mass is yet another factor which affects contact parameters. It was found 
that both stiffness and damping increased with increase in ball mass (see Table 
3.3). 

It has been reported by Venkataraman and Narayanan [11] that the impact 
force increases with ball diameter and is proportional to the square of the ball 
size. This explicitly points to the increase in stiffness with ball mass. Damping 
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Figure 3.6: Effect of particle layers on stiffness and damping values of materials 
for ball of mass 0.647 Kg dropped from a height of 0.25 m. 


Table 3.3: Effect of ball mass on contact parameters for drop height of 0.72 m and 
2 layers of magnetite particles. 



Ball mass (kg) 

0.05 

0.10 

0.25 

k 

4.7e7 

1.4e8 

2.52e8 

Q 

5e4 

17e4 

20e4 

&eq 

0.8e6 

2.8e6 

4.7e6 

Qeq 

180 

530 

720 
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being a function of mass and as well as a strong function of stiffness, also increases 
with ball mass (see Figure 3.7). 


3.6 Effect of Type of Materials to be Ground 


It was observed, for three types of materials, the stiffness and damping increased 
from quartz to limestone and further to magnetite (see Table 3.4). 

This trend can be attributed to structural flaws viz, pores and cracks. Materials 
having more flaws break easily i.e., they provide less damping, and have lower 
strength i.e., they possess lower stiffness (see Figure 3.8). 


Table 3.4: Effect of material type on impact parameters for ball of mass 0.647 Kg 
dropped on 2 layers of particles from a height of 0.25 m. 


Material type 

Quartz 

Limestone 

Magnetite 

k 

2.0e8 

2.4e8 

5e8 

<1 

1.3e5 

4.0e5 

5.0e5 

keq 

3.2e6 

4.0e6 

7.7e6 

Qeq 

509 

1249 

1897 


3.7 Determination of Coefficient of Restitution 

Coefficient of restitution is defined as a ratio between velocity just after impact, 
v a and velocity just before impact, v b . Mathematically, coefficient of restitution is 

. given as 


e = — — (3-1) 

Vb 

Calculation of coefficient of restitution e is of utmost importance because it is 
one of the input parameters required in simulation . In fact, a DEM code [18] is 
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Figure 3.7: Effect of ball mass on contact parameters for drop height of 0.72 m 
and 2 layers of magnetite particles. 





Section 3.7 Determination of Coefficient of Restitution 


40 




Material type 


Figure 3.8: Effect of material type on contact parameters for ball of mass 0.647 Kg 
dropped on 2 layers of particles from a height of 0.25 m. 
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available which takes material properties viz., stiffness, coefficient of restitution 
and coefficient of friction as input parameters. The values of parameters obtained 
by non-linear and equivalent linearized methods were used to calculate coefficient 
of restitution. The non-linear and linearized model were solved numerically using 
the corresponding non-linear and linearized stiffness and damping parameters. 
The velocity of ball just before impact was set at 2.4 m/ sec. Impact velocity of this 
magnitude is encountered in ball mills having diameter of 0.545 m where the ball is 
dropping from the mill centre vertically downwards on the mill wall. The impact 
process was assumed to be over when the impact force returned back to zero. Then 
the velocity just after impact was calculated. The ratio of final velocity to initial 
velocity with the sign changed gave the value of coefficient of restitution 1 . The 
following plots (see Figure 3.9 and Figure 3.10) show the variation of coefficient 
of restitution with damping for non-linear and equivalent linearized models. 



Non-linear damping, q 

Figure 3.9: Variation of coefficient of restitution with damping constant q for ball 
of mass 0.647 kg having impact velocity 2.4m/ sec (non-linear model). 

^ For linear model it is possible to arrive at the analytical solution for coefficient of restitution, 
for derivation see Appendix B. 




Cofficient of restitution. 
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Linearised damping, q 


Figure 3.10: Variation of coefficient of restitution e with damping constant q 
for ball of mass 0.647 kg having impact velocity 2 Am/ sec (equivalent linearized 
model). 




Chapter 4 

Simulation of a Ball Mill by DEM 


In this chapter a 54.5 cm diameter and 30.8 cm long ball mill is simulated by 
using a DEM based simulator [18]. Two types of contact models were used in the 
simulation viz . , linear and non-linear. However, the contact parameters in the case 
of a linear model were changed to incorporate the equivalent parameters extracted 
from the non-linear model. Several simulations were carried out and the computed 
power draw was compared with the corresponding experimental data. 

4.1 Mill Set-up for Simulation 

In this section, a laboratory size ball mill is considered to predict the power draw 
as a function of mill speed. This can be easily done by means of a DEM based 
simulator that tracks the motion of balls in all three directions. In order to as- 
sess the accuracy of the prediction power draw, available experimental data on a 
54.5 cm x 30.4 cm ball mill is considered. These data were published by Liddell 
and Moys [?], where the authors claim an accuracy of 0.5 percent in their data. 
Thus, these data are found most suitable for comparison purpose. 

A mill similar to that used by Liddell and Moys is considered for simulation. It 
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is a smooth mill without lifters. The grinding media is made up of graded charge 
of 0.029 m, 0.04 m, and 0.054m diameter balls having masses 0.096%, 0.252%, 
and 0.647 % respectively. To ensure mill filling of approximately 40%, 324 balls 
were generated randomly. In addition, the total load was maintained at 140 % by 
slightly increasing the density of the balls up to 8400 %/m 3 . 

The DEM based simulator [18], was fed with the above data and trial runs 
were carried out in SUN-Ultra450 machine. Figure 4.1 shows snapshots (front 
view) of the mill rotating at 60 % and 95 % of critical speed. The critical speed of 
mill is 6.05 radians /sec. 
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(a) 



(b) 


Figure 4.1: Figures showing snapshots (front view) of the mill for (a) 60% and 
(b) 95 % of critical speed. 
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4.2 Mill Simulation 

There are basically three contact stiffness models that are tested through the 
following numerical simulations. These are 

• Linear model 


F — kx 


(4.1) 


• Equivalent linearized model 

F = k eq x (4.2) 


• Non-linear model 


F = kx h6 (4-3) 

It is clear that all the three models are parameter dependent. In particular, the 
stiffness and the coefficient of restitution by and large decide the model response. 
Thus, at the very outset, the model parameters used in the respective simulations 
are shown in Table 4.1. This table shows the values of parameters that are used 
in the simulation such as stiffness, damping, coefficient of restitution, and coeffi- 
cient of friction. In the same table, the critical time-step required for the DEM 
simulation is also shown. These data are obtained as per the procedure described 
in Chapters 2 and 3. Briefly, the stiffness value in the case of a linear model was 
obtained by considering the initial slope of experimental force-deformation curve. 
In the case of non-linear model, the value of stiffness was directly extracted from 
the model (Equation 2.1) that matched the UFLC data. The equivalent param- 
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eters were then extracted from the non-linear model by applying the technique 
called equivalent linearisation. 


Table 4.1: Data used for simulation. 


Model type 

Linear 

Equivalent 

linearized 

Non-linear 

Normal stiffness k 


• 11.3e6 

5.53e8 

Shear stiffness k s 

3.0e5 

7.53e6 

3.69e8 

Coefficient of restitution e 

0.45 



Coefficient of friction /i 

0.90 

0.90 

0.2 

Critical time-step A t (sec) 

2.09e-4 

3.9e-5 



At first, the linear contact model was considered. The 0.545 m diameter mill 
was simulated by using the DEM based simulator where the contact model was 
appropriately set. At the end of the simulation, the power draw of the mill corre- 
sponding to the set of operating conditions was noted. To begin with, simulations 
were carried out for different values of coefficient of friction. The results obtained 
are shown in Table A. 5 in Appendix A. This was done in order to get the most 
accurate power draw for a single value of coefficient of friction. It was found that 
for a value of coefficient of restitution of 0.9, the maximum power draw predicted 
by the model was 93 % of the experimental maximum power draw. Finally, several 
simulations were conducted to establish the variation in power draw of the mill 
with its speed. All the simulations at each mill speed were carried out for three 
revolutions and an average power draw was considered for comparison purpose. 

Next, the non-linear and equivalent linearized models were tested in a similar 
manner. Since, the linear model predicted power draw close to 93 % for coefficient 
of friction 0.9, it was decided to use the same value of coefficient of friction for 
equivalent linearized model. Table 4.2 shows the power draw predicted by using 
the linear model, equivalent linearized model and the non-linear model for different 
mill speeds. For the purpose of comparison the experimental data is also presented 
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in the same table. The following observations are made from the table: 


• It is seen from the table that the non- linear model could not predict the 
experimentally observed power draw. 


Table 4.2: Comparison between simulated power draw results with the experimen- 
tal power draw. 


Mill speed 

Power draw (W) 

(% of critical) 

Experimental 

Linear model 

Equivalent 
linearized model 

Non-linear model 



H = 0.90 

H = 0.90 

H = 0.20 

50 

390 

384.2 

400.5 

547.5 

60 

490 

475.9 

488.6 

650.0 

70 

590 

546.0 

579.3 

753.7 

75 

640 

574.8 

627.7 

854.2 

80 

683 

610.4 

668.5 

1031.4 

90 

700 

651.1 

710.2 

1533.9 : 

95 

683 

630.6 

679.0 

968.8 j 


• The power draw predicted by equivalent linearized model agrees quite well 
with the experimental power draw. 

• Between the three models, the equivalent linearized model is found to be the 
most accurate one. 

The numerical data is also compared against the experimental data in Figure 4.2. 
It is seen from this figure that all three models show the same trend in the variation 
in power draw i.e., there exists a maxima in the power draw corresponding to a 
particular mill speed. Interestingly, all three model correctly predict this speed 
at which the power draw is maximum. However, the non-linear model could not 
match the power draw. It is believed that incorporation of the non-linear model 
into the 3-D DEM code has to be done more carefully. 



Power draw, 
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50 60 70 80 

Mill speed, % of critical 


Figure 4.2: Comparison between linear model, equivalent model and experimental 
results. 




Chapter 5 


Summary and Discussions 


The discrete element method (DEM) has become a powerful tool for systematic 
analysis of the charge behavior in many of the mineral processing units such as 
ball mill and jigs. However, it; suffers from the criticism that the discrete element 
models are highly parameter dependent and there is no unique way of establishing 
these parameters. This research work was undertaken to dispel this state of affairs 
by establishing a procedure to correctly determine the parameters viz., contact 
stiffness, damping, etc., while keeping the numerical and computational expense 
manageable. 

At first, the available experimental data on drop ball tests using a ultra fast 
load cell (UFLC) is considered. These data show the deformation behavior of the 
contact when a steel ball is dropped onto a steel anvil. • Analysis of these type of 
metal-metal impacts reveal that a non-linear contact model adequately represent 
the contact behavior. The contact model adopted was of the following form: 

mx + qx°' 8 x + kx 1 ' 6 = 0 (5.1) 


where, k is stiffness parameter and q is damping parameter. The parameters 
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involved in this model are extracted after matching the model response with the 
experimental data. These parameters allowed for the calculation of the energy 
involved during a given impact. This formed the basis for assessing the suitability 
of a particular contact model for the impact problem under investigation. 

The force-displacement behavior obtained during a collision (metal-metal) is 
modeled by using a non-linear as well as a linear contact model. It was found 
that the non-linear model best described the experimental data. In particular, 
the energy expended during the collision was very accurately predicted by using 
the non-linear model. However, the material stiffness as a parameter in the model 
turned out to be of the order of 10 10 N/m 1 ' 6 . This has serious ramification when 
this value is used in a discrete element model where the goal is to simulate the 
behavior of not just a single ball but thousands of balls. A large value of contact 
stiffness reduces the time-step of integration adding to the computational expense. 

In order to circumvent the above problem, a new approach is taken to transform 
the non-linear contact behavior to the most appropriate linear behavior in a formal 
sense. In this approach, the non-linear model is converted into its equivalent linear 
form using a equivalent linearisation technique. This method provided a linear 
stiffness value of 10 8 N/m which is two orders in magnitude lower than that of 
the corresponding non-linear stiffness parameter. Thus, the computational time 
significantly lowered by adopting the equivalent linearisation technique. 

In the foregoing, only metal-metal type collision problems are considered. But, 
it is realized that real ball mill simulation would involve particles that are to be 
ground. Under these circumstances, the contact model must incorporate material 
parameters that will take into account the presence of particles between colliding 
balls. Experimental data on drop ball test incorporating particle layers are gath- 
ered from the literature. Analysis of these data show that a contact stiffness of 
10 6 N/m sufficiently describes the contact behavior. This value of stiffness further 
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reduces the time-step in comparison to metal-metal collision problem. 

The performance of the equivalent linearized model is assessed by incorporating 
it into a generalized DEM software. A laboratory size mill of diameter 0.545 m 
and length 0.308 m was simulated. Three contact models of linear, equivalent 
linearized, and non-linear type were considered. Three types of balls were used 
in the simulation for which contact parameters were available. The ball mill was 
simulated to obtain the variation in power draw with mill speed. It was found 
while comparing with the available experimental data on a similar ball mill that 
equivalent linearized model predicted the power draw most accurately. 

The equivalent linearized model, although, is able to predict the power draw 
accurately, it fails to give the correct force-deformation and force-time histories. 
It is therefore suggested to use a displacement term in the relationship for the 
damping force so that initially when the deformation is zero the net force is also 
made equal to zero. 




Appendix A 


Experimental Data Used in the 
Present Study 


A.l Ultra Fast Load Cell Data Tables 

The tables below (Table A.l through Table A.4) depict the data collected by 
UFLC experiments which were used for parametric analyses. 


Table A.l: UFLC impact data with no particle layers in between. 


Experiment No. 

1 

2 

3 

Particle Layer 

0 

0 

0 

Ball Mass (kg) 

0.096 

0.647 

0.252 

Ball Height (m) 

0.30 

0.30 

0.30 

Impact Velocity (m/ sec) 

2.4 

2.4 

2.4 

Input Energy (J) 

0.28 

1.86 

0.73 

Consumed Energy (J) 

0.20 

1.08 

0.47 

Maximum Force ( N ) 

8109 

26000 

15025 

Maximum Displacement (mm) 

0.075 

0.136 

0.107 

Maximum Time (psec) 

78 

154 

116 
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Table A. 2: UFLC impact data on limestone. 


Maximum 
Displacement (mm) 

2.16 

Cl 

eo 

04 

OO 

04 

2.50 

co 

CO 

00 

ZVZ 

O 

Cl 

Cl 

cq 

CO 

r- 

CO 

xf 

Cl 

co 

04 

i q 
xF 

oo 

o 

in 

4.44 

3.16 

oo 

CO 

CO 

O 

o 

co 

rH 

o 

04 

1.74 

1.63 

O 

rH 

CO 

04 

xr 

CO 

xt 4 

04 

co 

Maximum 
Force (A'') 

1219 

2405 

7496 

17361 

2366 

04 

00 

xt 4 

00 

11277 

2587 

o 

CO 
xf 4 
CO 
1— 1 

370 

rH 

o 

CO 

m 

co 

o 

rH 


xT 

in 

co 

04 

CO 

CO 

in 



7783 


400 

1131 


Consumed 
Energy (J) 

0.92 

1.52 

2.95 

5.58 

1.52 

2.92 

5.89 

1.52 

5.83 

0.47 

0.98 




1 

1 

1 


j 


1 

1.17 

Input 

Energy (J) 

0.93 

1.55 

ire 

6.20 

1.55 

3.11 

6.20 

1.55 

6.20 

0.54 

1.08 


I 








I 

1.21 

Impact 

.Velocity (m/sec) 

1.70 

2.19 

3.10 

4.38 

2.19 

3.10 1 

4.38 

2.19 

4.38 

2.07 

2.92 

3.58 

4.13 

1.70 

2.19 

ore 

1.70 

2.19 

3.10 

1.70 

2.19 

3.10 

Ball 

Height (m) 

0.15 

0.25 

0.50 

O 

T— t 

0.25 

0.50 

O 

rH 

0.25 

O 

rH 

0.25 

0.50 

0.75 

1.0 

0.15 

0.25 

0.50 

0.15 

0.25 

0.50 

0.15 

0.25 

0.50 

Ball 

Mass (kg) 

0.647 

0.252 

0.647 

0.252 

Particle 

Layers 

CN 

LO 

O 

rH 

CO 

04 

m 

rH 

Impact 

Condition 

Ball/Flat 

Ball/Ball 
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Table A. 3: UFLC impact data on quartz. 


Maximum 
Displacement (mm) 

tH 

CD 

04 

00 

N- 

04 

2.10 











t-H 

F- 

o 

00 

CD 

Maximum 
Force (TV) 

04 

O 

00 

04 

04 

LO 

CD 

00 

CO 

04 

04 

OO 

t-H 

04 

CO 

04 

t-H 

LO 

r-H 

£>- 




CO 

o 

04 

t-H 

LO 

04 

O 

t-H 

O 

04 

O 

o- 

H 

04 

04 

F— 

F- 

04 

LO 

F- 

04 

04 

04 

CO 

04 

f- 

T i 

LO 

04 

F- 

CD 

04 

Consumed 
Energy (J) 

04 

t-H 

04 

CD 

04 

CO 

t-H 

LO 

LO 

04 

r-H 

04 

CO 

04 

CO 

t-H 

LO 

04 

O 

t-H 

CO 

LO 

t-H 

LO 

o 

04 

04 

t-H 

CD 

04 

LO 

F- 

F- 

04 

T— H 

CO 

XF 

04 

O 

oo 
















04 

LO 

LO 

Impact 

Velocity (m/sec) 

!>- 

CD 

04 

04 

04 

04 

CO 

r— i 

r- 

o 

04 

04 

04 

04 










CO 

t-H 

Ball 

Height (m) 

lO 

04 

O 

O 

LO 

o 

o 

o 

rH 

LO 

04 

CD 

LO 

O 

o 

q 

rH 

o 

LO 

CD 

LO 

o 

o 

o 

t-H 

LO 

04 

CD 

O 

LO 

o 

o 

o 

T— i 

LO 

04 

CD 

o 

lO 

CD 

o 

o 

rH 

Ball 

Mass (Kg) 

CO 

o 

04 

LO 

04 

O 

t- 

q 

o 

Particle 

Layers 

LO 

O 

T— 1 

CO 

LO 

O 

T 1 

Impact 

ondition 

E 

rH 

3 

"c3 

CQ 

ro 
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Table A. 4: ULFC impact data on magnetite. 


Maximum 
Displacement (mm) 

2.52 

05 

CQ 

csi 

r— ( 

05 

oa 

H 

CO 

o 

cvi 

00 

o 

csi 

LO 

0 

01 

rt 

rH 

Maximum 
Force ( N ) 

583 

1463 


CO 

CO 

977 

2715 

1145 

3709 

Consumed 
Energy ( J) 

0.29 

0.61 



1 



1.20 


I 





I 


1.24 

O 

co 

o If 

a . 

6 & 
t— i 0 

0 

1 

(N 

LO 

C5 

LO 

CO 

i> 

o 

LO 

oo 

LO 

H 

CN 

co 

t—4 

CO 

05 

CO 

T— 1 

CO 

C5 

i—l 

Ball 

Height (m) 

0.33 

0.67 

1.34 

0.13 

0.26 

0.51 

oro 

O 

<M 

O 

Ball 

Mass (kg) 

9600 

0.252 

0.647 

Particle 

Layer 

CN 

s 

-4-3 
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A. 2 Mill Simulation Data Tables 


Table A. 5: Power predicted by linear model for different values of coefficient of 
friction. 



Power draw (W) 

Mill speed 
% of critical 

Coefficient of friction 

0.3 

0.5 

0.7 

0.9 

60 

255.0 

332.7 

413.6 

475.9 

70 

298.9 

391.4 

480.7 

546.0 

75 

320.3 

418.5 

513.3 

574.8 

80 - 

346.2 

445.5 

534.8 

610.4 

90 

368.3 

467.6 

567.5 

651.1 

95 

354.8 

496.3 

557.1 

630.6 




Appendix B 


Important Derivations and 
Relations 

B.l Derivation of Damping Constant for Linear 
Model 

The equation of motion for the oscillating system consisting of a mass, spring and 
a dashpot is 


mx + qx + kx = 0 (B.l) 

The solution [16] of this equation under the initial condition that x = 0 and x = Vi 
at t = 0 is given by 

^ sin (rit) exp (-7 u 0 t) (B.2) 

'j {77 cos (rjt) - ju u 0 sin (r^t)} exp (—7 w 0 t) (B.3) 
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where 


\JkJm 

(B.4) 

q 

2 \frnk 

(B.5) 

“ 7 2 

(B.6) 


The oscillation period of this system is 2 ^/ 77 . A particle colliding with another 
particle at the time t — 0 detaches at time t = tt/t]. The velocity at the time 
t = -n/r] is given by 


Vo -i \t=,/ n = -Vi exp (7W 0 ir/i7) 


(B.7) 


Therefore, the coefficient of restitution becomes 


e = -v 0 /vi - exp (juoir/r)) 


(B.8) 


If the coefficient of restitution is regarded as a constant empirical parameter (it 
depends on impact velocity, in fact), the damping coefficient may be determined 
from the above equations. Consequently, the damping coefficient is given by 


Q = 



(B.9) 


where m is the ball mass, k is the contact stiffness, t is he impact time, and x is 


the ball deformation. 



Section B.2 Relation Between Coefficient of Restitution and Impact Velocity 


60 


B.2 Relation Between Coefficient of Restitution 
and Impact Velocity 

It was observed by Goldsmith [9] that the magnitude of coefficient of restitution 
decreased monotonously from unity with increasing impact velocity (Figure B.l). 
Also, the reduction in value of coefficient of restitution was greater than can be 
justified by the loss of energy due to damping. This feature contradicts the basic 
assumption of the Hertz theory of contact regarding perfect elasticity of process 
in which in absence of damping force the coefficient of restitution becomes unity. 
Therefore, it is revealed that apart from energy loss due to damping, there is some 
extra loss of energy due to plastic deformation of balls. It follows that the model 



Figure B.l: Coefficient of restitution as a function of impact velocity for spheres 
of same size but different materials, after Goldsmith [9]. 

developed in the present work requires some modifications to take into account 
the permanent deformation of balls by taking into consideration the strain-rate 
effects in the material behavior. 
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